preprint 



Diffusion rates of Cu adatoms on Cu(lll) in the presence of an adisland nucleated at 

FCC or HCP sites 

Mihai-Cosmin Marinica*, Cyrille Barreteau*, Daniel Spanjaard^ and Marie-Catherine Desjonqueres* 
"CEA Saclay, DSM/DRECAM/SPCSI, Bailment 462, F-91191 Gif sur Yvette, France and 
' Laboratoire de Physique des Solides, Universite Paris Sud, Batiraent 510, F-91405 Orsay, France 

(Dated: February 2, 2008) 

The surface diffusion of Cu adatoms in the presence of an adisland at FCC or HCP sites on 
Cu(lll) is studied using the EAM potential derived by Mishin et al. [Phys. Rev. B 63 224106 
(2001)]. The diffusion rates along straight (with close-packed edges) steps with (100) and (lll)-type 
microfacets (resp. step A and step B) are first investigated using the transition state theory in the 
harmonic approximation. It is found that the classical limit beyond which the diffusion rates follow 
an Arrhenius law is reached above the Debye temperature. The Vineyard attempt frequencies and 
the (static) energy barriers are reported. Then a comparison is made with the results of more 
realistic classical molecular dynamic simulations which also exhibit an Arrhenius- like behavior. It 
is concluded that the corresponding energy barriers are completely consistent with the static ones 
within the statistical errors and that the diffusion barrier along step B is significantly larger than 
along step A. In contrast the prefactors are very different from the Vineyard frequencies. They 
increase with the static energy barrier in agreement with the Meyer-Neldel compensation rule and 
this increase is well approximated by the law proposed by Boisvert et al. [Phys. Rev. Lett. 75 
469 (1995)]. As a consequence, the remaining part of this work is devoted to the determination of 
static energy barriers for a large number of diffusion events that can occur in the presence of an 
adisland. In particular, it is found that the corner crossing diffusion process for triangular adislands 
is markedly different for the two types of borders (A or B). From this set of results the diffusion 
rates of the most important atomic displacements can be predicted and used as input in Kinetic 
Monte-Carlo simulations. 

PACS numbers: 68.43.Fg, 68.43.Hn, 68.43.Jk 



I. INTRODUCTION 

The diffusion of adatoms on metal surfaces is still the subject of very active research 

HQ 

. Indeed, it plays a crucial 

role in crystal growth which is important to master in view of the applications, for instance in nanotechnologies. The 
more the growth is understood at an atomic level, the more the fabrication processes can be controlled in order to 
obtain a better quality of the device performances. 

As a starting point, homo-epitaxial systems can be used as models since effects such as lattice mismatch or inter- 
diffusion processes are excluded. A large number of theoretical and experimental works has been devoted to homo- 
epitaxial growth, e.g., on Pt (ESQ HQ, Cu (sUlQ, A S QQ, Rh QQQ, Al QQQ In this respect, the 



understanding of adisland shapes is of fundamental importance since they can indirectly influence the growth mode. 
The (111) surface of FCC metals is particularly interesting due to the observation of fractal as well as two (2D) or three 
(3D) dimensional compact adislands depending on the temperature [sUfiillfj. Actually, the adisland morphology and 
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its evolution with temperature result from a competition between thermodynamics and kinetic phenomena depending 
on the activation (or not) of various diffusion processes. The shape of adislands is governed by thermodynamics when 
the temperature is such that the adisland is able to relax to its equilibrium configuration during the time interval 
separating the incorporation of two consecutive atoms. In this case the shape of adislands on a (111) FCC surface is 
determined by the ratio of two step energies since two types of steps with close-packed edges exist on this surface: 
the A step with a ledge of (lOO)-type orientation and B step with a ledge of (lll)-type. The anisotropy ratio of the 
two step energies / -E^gp is most often very close to unity (for instance, using the Mishin et aL[2j| potential for 

Cu, we have found i?^. =263meV and £^ e p=265meV [21]]) and therefore the adisland shape derived from the Wulff 
theorem is an almost regular hexagon exhibiting three-fold symmetry with edges corresponding to the smallest step 
energy slightly longer than the other ones [3- When kinetic effects dominate, the adisland morphology depends on 
the flux of atoms impinging on the surface as well as on the rates of adatom diffusion on terraces, along steps, around 
adisland corners and of adatom attachment to (detachment from) the adisland. 

The most efficient tool to simulate the evolution of the adisland morphology is the Kinetic Monte-Carlo (KMC 
[31313) method which needs as input the various rates of all the elementary atomic processes. Two types of 
processes are essential in this respect: the diffusion along steps (and its anisotropy between step A and step B), and 
the corner crossings to go from one edge to another. At low temperature when none of these processes are activated 
one expects a fractal growth if the impinging flux is not too small since diffusion is rapid on a flat (111) FCC surface 
(e.g. the activation energy is only 40meV for Cu). When these two processes occur, compact adislands can grow. 
An intermediate situation is also possible corresponding to an activation of step diffusion while corner crossing is 
still frozen. However other physical properties may have some importance: for instance, the step on which adatoms 
or dimers preferentially bind and the kink dissociation process play a role on the adisland shape which can vary 
with temperature from triangles limited by A steps, to triangles limited by B steps with intermediate hexagon-like 
patterns 0, 0] . The various possible sequences are extremely dependent on the details of the energy profile and the 
modification of a single barrier can lead to very different adisland shapes. Moreover the diffusion across a step plays a 
crucial role in the crystal growth mode (Volmer- Weber, Stransky-Krastanov or Frank Van der Merwe). Indeed during 
growth adatoms can be deposited either on a wide terrace or on a preformed adisland. In the former case, the atom 
will diffuse until it sticks to the adisland, whereas in the latter case the adatom needs very often to overcome an 
extra activation energy (Schwoebel barrier) to incorporate to the descending step edge. Depending on temperature 
a 3D growth can thus be initiated. Therefore we have carried out a very detailed analysis of the various possible 
diffusion processes of a copper adatom in the vicinity of an adisland. Two types of epitaxy will be considered in which 
adatoms occupy either normal FCC sites or faulted HCP sites. Actually the stacking fault energy is very small in 
copper and adisland nucleation has been observed at FCC sites on the nominally flat surface and at HCP sites on the 
Cu(21, 21,23) vicinal surface 

Many experiments based on field ion|3 13 anc ^ scanning tunnelling 0, |3 13 microscopies have been devoted 
to the determination of surface diffusion coefficients. However their interpretation may be somewhat tricky since, 
apart from peculiar cases, several elementary processes are involved. On the theoretical side, diffusion rates have 
been deduced either from Transition State Theory (TST) in the harmonic approximation |2lJ, l3fi l3ll |32j (TST-HA) , 
or using classical thermodynamical integration |3J| (TI) (see Appendix A) or else from classical molecular dynamic 
(MD) simulations |34l l35l l36l 1371 138| . In the classical limit all theoretical approaches conclude that an Arrhenius law, 
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r = To cxp(— AE/ksT), fits accurately the evolution of the diffusion rates with temperature in good agreement with 
experiments. On the basis of MD simulations it has been suggested that the barrier AE for diffusion along steps 
on Ag(lll) and Au(lll) might be different from the static barrier whereas Boisvert et al. found that the 
value of the static barrier for Cu on Cu(100) lies inside the error bars of the results obtained both with MD and 
TI methods. Actually a very good accuracy on AE (better than a few 10 _2 eV) would need huge simulation times. 
Furthermore the Arrhcnius law being obviously approximate, the separation between a prefactor and an exponential 
term is somewhat arbitrary as will be discussed in Scc.IIB. 

From the above remarks, it is clear that any investigation of diffusion processes should begin with the determination 
of the (static) potential energy barriers. A limited number of barriers concerning high symmetry systems have been 
calculated using ab-initio codes but these codes are too computer time demanding when the symmetry is low, 
i.e., when the number of elementary diffusion processes is large. One must then rely on semi-empirical potentials 
such as those derived from Effective Medium Theory (EMT)j^|, Embedded Atom Model (EAMl [4li or Second 
Moment Approximation (SMA)^J. Even with these simple models, unknown barriers are sometimes deduced using 
some approximations Q . Even though these approximations are reasonable some effects can be missed as we will see 
in the following. 

The aim of this work is the determination of the diffusion rates corresponding to most of the elementary diffusion 
processes that may occur during the growth of a single 2D adisland on Cu(lll). After having briefly presented the 
potential used (Sec. HA) we study the diffusion of an adatom along step A and step B (Sec. Ill) using the two methods 
described in Sec.IIB, namely the TST-HA approach (in the classical limit) and MD simulations. We show that the 
static activation barrier accounts quite well for the results of MD simulations. However, the prefactors obtained in the 
latter method are different from those given by TST-HA. Furthermore they are quite consistent with the Meyer-Neldcl 
law proposed by Boisvert et al. js^ . As a consequence, in the remaining part of the paper, we limit ourselves to 
the determination of the static barriers for other diffusion events in the presence of straight steps (Sec. IV) or in the 
vicinity of steps with defects, for instance around corners and kinks (Sec.V). Conclusions are drawn in Sec. VI. Finally 
Appendix A briefly summarizes the TST-HA theory and the TI method. 

II. FORMALISM 
A. The potential 

I — I 

We have used in this work the EAM potential derived by Mishin et aL|2£| for copper. In this model the total 
energy of an assembly of N atoms with respect to that of N isolated atoms is written as a function of all interatomic 
distances r^: 

N N 
i,j=l i=l 
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where pi can be interpreted as a function proportional to the electron density induced at site i by the neighbors, i.e.: 

Pi = ^2p(nj)- (2) 

i¥* 

The proportionality factor is chosen such that pi — 1 for a bulk atom at equilibrium. The chosen reference energy 
implies that V(r) and p(r) vanish when r tends to infinity and F(0) = 0. The functions V(r), p(r), and F(p) arc fitted 
parameterized functions which, on the whole, contain 26 parameters. The 26 parameters are required to give exactly 
the bulk equilibrium lattice parameter, cohesive energy and bulk modulus and to fit selected properties of copper 
taken from experiments or obtained by ab-initio calculations. These properties are mainly related to the bulk phase 
and include in particular the stacking fault energy which is of a peculiar importance in the present work. More details 
can be found in Ref. In addition, Mishin et al. have shown that their potential was able to reproduce many 

physical quantities not included in the fitting data base and, in particular, the surface energies of low index surfaces. 
Moreover we have already used this potential in a previous work [21I ]. First we verified that it gives other surface 
properties, like the step and kink energies, with numerical values very close to experimental data. Then we used it to 
study the diffusion of monomers, dimers and trimers on Cu(lll) and our results were in very good agreement with 
STM observations. We are thus very confident in using this potential for the present problem. 

B. Determination of diffusion coefficients along straight steps 

The diffusion coefficients along straight steps (i.e., with a close-packed edge) have been calculated using two tech- 
niques. On the one hand, they can be derived from the TStJ^J] by combining the determination of the minimum 
energy path and the calculation of the vibrational free energy in the framework of the HA and, on the other hand, 
using MD simulations. Indeed, as will be seen in the following, the most frequent diffusion event for an adatom 
attached to a straight step corresponds to jumps along this step since the barriers corresponding to other events are 
much higher. Thus this diffusion is essentially a one-dimensional problem. 

1. Diffusion coefficients of an adatom along straight steps from TST-HA 

From TST the surface diffusion coefficient D(T) along straight steps at temperature T is related to the diffusion 
frequency (or diffusion rate) Ttst(T) by: 

D(T) = ia 2 r TST (T) (3) 

when the diffusion is assumed to proceed by uncorrelated jumps of length a between neighboring adsorption sites 
along the step (note that a increases slowly with T when thermal expansion is taken into account but in practice this 
effect is quite negligible). The diffusion frequency Txst(T) is given by: 



r TST (T) = 2^ exp(-AF/fc B T) (4) 
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where AF is the difference of total free energy between the saddle point and the starting adsorption configuration 
(stable or metastable equilibrium). The factor 2 arises from the number (n c = 2) of diffusion channels, i.e., Ttst(T) 
is the diffusion frequency in both directions along the step. AF can be split into two contributions: 

AF = AE + AF mb (5) 

AE is the potential energy barrier, i.e., the difference in total energy (EqQ between the saddle point and the 
equilibrium configuration. This static diffusion barrier AE can be calculated by determining the minimum energy 
path using the Ulitsky-Elber alg orithm Q] which we already used with success in our study of the surface diffusion 
on the flat Cu(lll) surface |21|. For technical details the reader is referred to this latter reference. The quantity 
AF V ib = AUyib — TAS V ib is the contribution of vibrations to the variation of the free energy (with obvious notations). 
In the HA (see Appendix), AF v n, is given by: 

AF vib = k B T ln(2sinh(— — ))An(i/)<^ (6) 
Jo 2fc B T 

where An(v) is the difference between the vibrational densities at the saddle point and at equilibrium which are 
computed from the eigenvalues of the dynamical matrix. Note that at the saddle point one of the eigenvalues is 
negative. This eigenvalue is excluded from the vibration density and, consequently: 



An{v)dv = -1. (7) 



o 

Thus in this TST-HA model, the diffusion rate is: 



u,.e 



r?s\(T) = t^-S e xp(-AE/k B T) = Vg A (T) e x p(-AE/fc fl T) (8) 



3JV-1 , s 

hv s p 



n ^2k B T' 
p=i 

if we call and i/| the (real) eigenfrequencies for the equilibrium and saddle point configurations and N + 1 the total 
number of atoms, i.e., including the adatom. Note that the three free translational modes (of zero frequency) have 
been excluded from the spectra. It is seen that the prefactor Tq A (T) depends in general on T and, consequently, 
^tst(T) does not follow an Arrhenius law. 

Let us examine the limits of low and high temperatures. In the low temperature limit AS v n, —> and AU v n,(T) = 
AU mb (0) + 0(T 4 ) Q so that: 



r 



k R T 



£# T (T) ~ 2^|i exp(-(AE + AU mb (0))/k B T) (9) 



In the high temperature limit we get: 

r Tsr( T ) = 2iy o cxp (-AE/ksT) (10) 
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in which v$ is the Vineyard attempt frequency given by |47| : 

3N 

n 

= 5^5— (11) 

p=i 

Thus, in this limit, T^g T (T) follows an Arrhcnius law with a diffusion barrier AE and a prefactor 2z^o- 

In this TST-HA approach the deviation from the Arrhenius law comes from quantum effects. We will see in the 
following that these effects become negligible above the Debye temperature (To = 343 A' for Cu Accordingly the 

classical molecular dynamics approach, which we will briefly describe in the next subsection, is fully justified above 
T D . 



2. Diffusion coefficients along straight steps from MD simulations 

The motion of the N+l atoms is studied by MD simulations in which the Newton equations of motion in the 
potential given by Eq. ^ is solved using the Verlet algorithm in its velocity form. The simulation is carried out 
in the temperature range 350K-600K and no correction of the temperature is needed since quantum effects can be 
neglected [49j . The thermal expansion of the Cu crystal is derived from MD simulations in the bulk and the nearest 
neighbor distance is fitted by a = 2.553 + 2.51610~ 5 T + 1.193910~ 8 T 2 in A. After an equilibration period, the motion 
of the adatom is observed during another period t run and the trajectory of the adatom is recorded. This trajectory 
is analyzed by assuming that the diffusion proceeds by uncorrelated and discrete jumps. At this point, we must 
emphasize the differences between MD and TST. The basic assumption in TST is that each crossing of a divi ding 
surface containing the saddle point corresponds to a diffusion event and thus "recrossing effects" are ignored |44j . 
This is not the case in MD since recrossing events can be identified, i.e., they occur when the adatom turns back in 
a very short time after having crossed the dividing surface. Moreover in Eq|3]each diffusion event is assumed to be a 
jump of length a while long jumps of length na can occur in MD. Finally, anharmonic effects are included since the 
forces are calculated from the exact expression of the potential. 

The diffusion coefficient is given by the Einstein formula: 

Dm _ to <«Efl> ,12, 

where 

N S 

r(t) = ^riiYi, (13) 

Ns is the total number of jumps during the time interval t, riiTi is the displacement vector corresponding to the i th 
jump (rij = 1,2,3... for a simple, double, triple... jump, respectively, and = a) and, for uncorrelated jumps (i.e., 
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< ivrj >= 0), if t run is large enough: 



D{T) 



a z vA n 



(T) 



(14) 



where N n is the number of jumps of length na observed during the period t Tun - 

The discrimination between the accepted and rejected diffusion events and between the different jump lengths is 
done with the same criterion as in the work of Ferrando and Treglia js^ . Let us call r t h the thermal time given by: 



Tth 



1/2 



(15) 



where m is the mass of the adatom, i.e., r t ^ is simply the time to go over the lattice spacing for a particle with an 
energy fc^T, for example r t h = lps at T=500K. The diffusing atom is considered to have resided at a given site if 
the time spent on this site is longer than Tth- This is illustrated in FigQ For instance the case (a) corresponds to a 
single jump since the time spent at sites i — 1 and i is longer than T t h, whereas no jump is recorded in the case (b) 
since the time spent at i + 1 is too short (recrossing event). Finally the case (c) corresponds to a double jump since 
the time spent at site i is too short. 
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FIG. 1: Typical time evolution of the position of an adatom along a step. The adatom starts from site i and the unit along the 
ordinate axis is the nearest neighbor spacing a. The adatom is considered to be at site i if its coordinate lies in the interval 
[i — 1/2, i + 1/2]. The rejected, single and double jumps are marked by arrows. 

The usual method to analyze the temperature dependence of Y MD is to draw an Arrhenius plot, i.e., \nT MD {T) 
vs l/ksT jsil I34I I35I Isfii I37I Issj. Most often a quasidinear behavior is observed, the larger the number of recorded 
diffusion events (i.e., the smaller the statistical error), the smaller the deviation from linearity. As a consequence, 
yMD (rp-^ can a pp roacncc i by an Arrhenius law: 

r™(T) = r™ cxp (-AE MD /k B T) (16) 
in which the parameters Tq ID and AE MD are obtained by a least mean square fit. 
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It is then found that AE MD differs slightly from the static barrier AE. This small discrepancy may be due 
to statistical errors. Actually, as we will illustrate in the following (see FigEJl, if the fit is carried out by setting 
AE = AE the corresponding straight line is also contained inside the error bars. Furthermore, if real physical 
effects were responsible for the variation between AE MD and AE, then (AE MD — AE) should be a function of 
temperature. Consequently the Arrhenius law is not strictly obeyed and the splitting of r MD into a prefactor and 
an exponential is somewhat arbitrary. However if the variation of (AE AID — AE) is of the first order in T, then the 
Arrhenius law remains almost strictly obeyed with a slope equal to AE. 

Actually, using TI (see Appendix) Boisvert et al [3^ have shown that the diffusion frequency can be written as 
an Arrhenius-type law in which the prefactor Tq 1 and the barrier AE TI are almost independent of temperature in 
the range 100 — 800K for surface diffusion of Cu on Cu(100). Moreover, by comparing the TI method with MD 
simulations, these authors observed that the static barrier AE lies always inside the error bars of AE TI and AE MD . 
We will see in the following that, similarly, our results on the diffusion of Cu along straight steps of Cu(lll) can be 
fitted nicely by assuming AE AID = AE. 



III. DIFFUSION OF A CU ADATOM ALONG STRAIGHT STEPS ON CU(lll) 



A. The geometry 

It is well known that there are two types of steps with close-packed edges on Cu(lll): the step A with a (100) ledge 
and the step B presenting a (111) ledge. The study of the surface diffusion along or across these steps is a prerequisite 
to investigate the growth of large Cu adislands on Cu(lll). These adislands are expected to be bordered by steps of 
type A or B and, out of equilibrium, may be or not in stacking fault|lo| since the stacking fault energy is rather small 
in copper. In the first case the atoms occupy HCP sites while, in the second, they are located at FCC sites. Thus we 
have to study four geometries: step A or step B for adislands at FCC or HCP sites. 
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FIG. 2: Diffusion processes of a Cu adatom (gray circles) along the steps A and B limiting a stripe in FCC geometry on 
Cu(lll): j (jump) and x (exchange) mechanisms. The filled black circles are the atoms of the stripe and the atoms of the 
substrate are denoted as open circles with a size decreasing when going towards the inside of the substrate. 
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In practice a super-cell containing pi x p2 two dimensional (2D) unit cells with N p layers is built to represent the 
Cu(lll) surface. On this super-cell a stripe of N s t r i P e close-packed atomic rows is added at FCC or HCP sites. This 
stripe is thus limited by an A step on one side and a B step on the other one (Fig[3Jl. The number N str i P e of rows is 
large enough to avoid any interaction between the two steps. On each side of this stripe (A and B edges) an adatom 
is deposited. Finally, usual periodic boundary conditions are applied. 

B. Diffusion of adatoms in the TST-HA approach 

We have first carried out static calculations from which the minimum energy path between two first neighbor 
adsorption sites along A and B-steps is determined for the jump and exchange diffusion mechanisms shown in Fig|2 
These calculations were performed on a large super-cell: the slab representing the surface contained (11 x 11) 2D unit 
cells and was 10 layer thick. The stripe was made of 4 atomic rows. In view of the large number of atoms in this 
super-cell, only the point k// = in the surface Brillouin zone was used to calculate the phonon frequencies. Similarly 
to our previous work the equilibrium structure was deduced from the conjugate gradient method and convergence 
was achieved when all the forces were smaller than 10~ 3 eV/A. The minimum energy path was determined using the 
Ulitsky-Elber algorithm [4^ and the iteration process was stopped when all the forces perpendicular to the path were 
less than 31CT 2 eV/A. 



Step type 


A 


B 


Stripe geometry 


FCC 


HCP 


FCC 


HCP 


AEj (meV) 
T$ A = 2i/ (THz) 
AE e (meV) 


247 
7.84 
1532 


235 
8.54 
1512 


312 
7.44 
1715 


302 
8.04 
1691 



TABLE I: Diffusion barrier AEj and corresponding prefactor (vq is the Vineyard attempt frequency) for Cu self-diffusion 
(jump mechanism) along A and B-steps on the (111) surface, in regular FCC and faulted HCP geometries. In the last line we 
have also reported the diffusion barrier AE e for the exchange mechanism shown in Fig|5| 

The results show that the adsorption energy of an adatom along a B-step is favored by 7 meV with respect to 
an A-stcp. This difference in energy is probably too small to have any influence on the island growth morphology 
and could not be at the origin of an asymmetry between A and B-steps. In contrast the diffusion paths along step 
A and step B for the jump mechanism are quite different. While along step A the diffusion proceeds by FCC — > 
HCP(saddle point)— > FCC jumps on the (111) terrace, along step B the path can be approximated by a displacement 
along the channel of a (Oil) microfacct which can be considered as the ledge. Indeed let us recall that the (lll)FCC 
surface with periodic monoatomic B steps can be denoted either as (p+l)(lll) x (111) or as p(lll)x(011)^|. The 
resulting energy barriers presented in Table[I]show that diffusion always proceeds by the jump mechanism (the energy 
barrier for the exchange mechanism being extremely high) and is expected to be faster along step A than along step 
B since its energy barrier is smaller by about 60mcV. Moreover, for both steps the diffusion barrier along a HCP 
stripe is slightly smaller by lOmcV than along a FCC stripe, i.e., the variation of the energy barrier between A and B 
steps is almost exactly the same in FCC and HCP geometries: 65meV and 67meV, respectively. Consequently from 
static calculations, the asymmetry between the diffusion along A and B-steps should be rather similar in FCC and 
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HCP geometries. Then we have calculated the phonon spectrum for adatoms at the equilibrium and saddle point 
configurations, from which the prefactor Tq A (T) is deduced. In FigOt we show the evolution with temperature of this 
prefactor: it is seen that the high temperature limit is attained above ~ 300K, i.e., close to the Debye temperature of 
Cu. This limit is equal to twice the Vineyard attempt frequency fEa lll|) and is given in Tabled Accordingly, above 
300K, the diffusion rate follows an almost perfect Arrhenius law which is actually almost undistinguishable from its 
asymptotic limit 2vq cxp(— A_E/fcsT) (Fig|3)). One notes that the prcfactors are not very dependent on the geometry 
and thus, in this model, the diffusion rate is governed almost entirely by the exponential term. The diffusion along 
step A is indeed faster than along step B and the motion is slightly easier in the HCP geometry. In order to quantify 
the asymmetry between step A and step B we have also calculated the ratio Ya/^b of the diffusion rates along the 
two types of steps (see Fig. 0J. This ratio is very similar in HCP and FCC geometries, it increases as the temperature 
decreases and typically adatom diffusion processes are approximately 4 and 6 times faster along step A than along 
step B at 600K and 450K, respectively. 
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FIG. 3: a) Prefactor Tq (T) of the jump rate, obtained from the TST-HA model, for Cu diffusion along A and B-steps on 
Cu(lll), in FCC and HCP geometries, b) Arrhenius plot of V^g T (T) corresponding to the four geometries in the temperature 
range 340K - 645K. 



C. Diffusion of adatoms from MD simulations 

The simulations have been performed on a (11 x 11) slab containing 16 layers with a stripe made of 4 atomic rows in 
FCC sites on one side of the slab and in HCP sites on the other side. One adatom is deposited along each of the four 
close-packed step edges of the two opposite stripes. This procedure allows to obtain with a single simulation box, the 
four trajectories from which the diffusion coefficients are extracted. We have performed a scries of MD simulations 
for temperatures ranging from 350K to 600K. The classical equations of motion were integrated with a time step of 
3.5fs. The system was first equilibrated for about 12ps. In the considered temperature range the far most frequent 
event is a jump along the step, however from time to time an extra event occurs, such as the formation of a dimer 
with an atom escaping from the step (see Sec. IVBjl . To avoid the occurrence of such "unwanted events", we have 
preferred to perform a number of independent simulations rather than a single and longer one. At each temperature 
the total simulation time was chosen in order to get a meaningful statistics, i.e., it increases when the temperature 
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FIG. 4: Arrhenius plots of the ratio Fa/Tb of the jump diffusion rates for a Cu adatom along A and B steps on Cu(lll) 
in FCC and HCP geometries as deduced from the TST-HA model and MD simulations (open circles: FCC geometry, open 
triangles: HCP geometry). 

increases. The diffusion coefficient can be calculated either from the Einstein relation (Ea ll2[) . or from Eall4l which 
assumes uncorrclatcd discrete jumps. We have found that the two methods give practically the same results showing 
that the jumps are indeed uncorrelated. In the following the results will be analyzed using the second approach which 
has the advantage of allowing an estimation of statistical errors. Thus the analysis of the trajectory consists in an 
enumeration of the various jumps (accepted or rejected) occurring during the simulation time. The statistics of the 
number of jumps is presented in Table [HJ It is seen at first glance that significant differences with TST-HA are 
expected since the transmission coefficient k (i.e., the fraction Ni/(N± + N r ) of accepted jumps) is far from being 
unity. 



Step type 


A 


B 


Stripe geometry 


FCC 


HCP 


FCC 


HCP 


T(K) 


t(ns) 


Ni N 2 N 3 N r 


Ni N 2 Nz N r 


Ni N 2 N 3 N r 


Ni N 2 N z N r 


350 
400 
450 
500 
550 
600 


165.9 
98.0 
77.7 
58.8 
39.2 
30.8 


277 3 258 
391 17 1 382 
696 21 1 582 
960 33 3 875 
977 55 4 993 
1159 57 7 1020 


401 9 180 
608 19 1 343 
953 30 2 553 
1294 38 3 727 
1323 60 3 815 
1448 88 9 926 


71 137 
127 3 1 259 
275 8 478 
447 19 779 
540 26 4 842 
709 38 7 1249 


115 3 82 
193 7 187 
377 16 1 384 
544 31 5 561 
603 50 4 722 
769 57 8 910 



TABLE II: Statistics of single (N\), double (N 2 ), triple (Nz) and rejected single (N r ) jumps for the diffusion of a Cu adatom 
along step A and step B on Cu(lll), in FCC and HCP geometries, from MD simulations. The total simulation time t is given 
for each temperature T. 

In Fig|S] we show the Arrhenius plots of T MD (T) for the four geometries with error bars corresponding to the 
standard deviation. We have first carried out a least mean square fit with the two parameters Tq ID and AE MD . 
Similarly to the results of Boisvert et al. 33] on the surface diffusion of Cu on Cu(lOO), we find that AE AID is quite 
close to the static barrier AE (see Table llllf) and the fitted straight line lies inside the error bars demonstrating the 
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quality of the statistics. However, as already discussed in Sec.IIB2, it is advisable to compare with a fit in which the 
barrier has been fixed to its static value. The corresponding fits shown in Fig. 0are nearly as good as the previous 
ones. Moreover the associated prefactors become quite comparable in FCC and HCP geometries in contrast with the 
previous fits for which the prefactor for the B step in the HCP geometry was much smaller than in the FCC geometry. 
This does not seem physically reasonable but is most probably due to the uncertainty on AE MD since a small error 
in AE MD must be compensated by a large variation of the prefactor in order to give the same value of T MD . Thus 
we have adopted the values of the prefactors Tq ID derived from the fit in which AE MD = AE. In this case it is 
seen from Tables HI and IIIII that Ta/^b as deduced from MD is about twice smaller than the corresponding value in 
TST-HA (see FigQJ. Indeed in the TST-HA model the prefactors are very similar for steps A and B while, from MD 
simulations, the prefactor corresponding to steps B, which have the highest diffusion barrier, is about twice that of 
step A. We will now give a physical interpretation for this difference. 



Step type 


A 


B 


Stripe geometry 


FCC 


HCP 


FCC 


HCP 


< K > 


0.52 


0.64 


0.36 


0.50 


(a) 


AE MD (meV) fitted 


239 


226 


308 


279 




r£ /D (THz) fitted 


5.0 


4.9 


11.6 


7.8 


(b) 


AE (meV) fixed 


247 


235 


312 


302 




r£ /D (THz) fitted 


6.0 


6.1 


12.8 


14.0 



TABLE HI: Transmission coefficients < k > averaged over temperature, energy barriers and prefactors derived from MD 
simulations: a)both parameters AE MD and T^f D have been obtained from a least mean square fit; h)YQ ID has been obtained 
from a least mean square fit with AE MD fixed at its static value AE. 



It is frequently found that when the activation energy increases within a family of processes, the prefactor also 
increases. Thus, in the expression of the process rate, the increase of the prefactor somewhat "compensates" for the 
decrease in the Arrhenius exponential term governing the dependence on temperature. This effect is known as the 
Meyer- Neldel compensation law jjn] . It has been investigated by Boisvert et al. 35] in the particular case of surface 
self-diffusion by means of MD simulations. By studying jump and exchange processes on the (100) and (111) surfaces 
of elements belonging to the end of the transition series or to noble metals, they have proposed the following law for 
the prefactors: 



r = r 00 exp(AE/A ) a (17) 

From theoretical models the exponent a is expected to lie in the range 0.5-1 depending on the nature of the 

excitations that give rise to the activated process. For acoustical phonons, a simple phcnomenological model predicts 
a — 3/4 and a characteristic energy Ao of the order of 2-3 times larger than a typical phonon energy (for instance 
hu = ksT]j where To is the Debye temperature). Actually Boisvert et al. have found that their results can be 
reasonably fitted with a = 0.7, A = 74meV and r 00 = 0.74THz. 

As already stated, our MD data for steps A and B obey the same trend since T^f D increases with AE. Thus it is 
interesting to add the corresponding data (see Fig(HJ) on the plot giving the logarithm of the prefactor versus (AE) 7 
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FIG. 5: Arrhenius plots of the MD jump rates for Cu diffusion along straight steps on Cu(lll) in FCC and HCP geometries. 
The error bars give the statistical error. The straight line is a least mean square fit in which the energy barrier has been fixed 
to its static value, a) Step A in FCC geometry; b)Step A in HCP geometry; c)Step B in FCC geometry; d)Step B in HCP 
geometry. 

which should be a straight line (Fig. 3 of RefQ. Our results are quite close to the straight line obtained by Boisvert 
et ai. 
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FIG. 6: Plot of In To as a function of AS ' 7 . The points represented by empty triangles are taken from Boisvert et al. |35l and 
refer to various surface self diffusion processes on low index surfaces of Au, Ag, Ni and Pd. The points indicated by arrows 
(filled black diamonds) are derived from our molecular dynamic simulations along A and B straight steps on Cu(lll) in FCC 
and HCP geometries. The dashed line is the fit corresponding to Eq |17l 
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Let us now comment on the differences between the TST-HA and MD approaches. As already mentioned, MD 
simulations take into account the transmission coefficient k and anharmonic effects. From our simulations k does not 
vary significantly with temperature. From Table ITTT1 it is seen that on average ka/kb — 1.5 for both geometries. As 
a consequence the increase of the prefactor for step B with respect to step A is not due to the transmission coefficient 
but rather to anharmonic effects which, as discussed by Boisvertei al. |.S j. Isjjj . lead to multiphononic excitations. The 
prefactor is then proportional to the number of ways of assembling these excitations. This gives rise to an entropy 
factor increasing with the barrier height (sec Appendix). This interpretation is consistent with the absence of a 
Meyer- Neldel effect in the TST-HA model which has already been noted by other authors j^l ■ 

As a conclusion, it is clear that the prefactor should increase with the (static) barrier height. Moreover the law 
proposed by Boisvert et al. (Ea ll7| can be used to get a reasonable estimate of this variation. Thus, in the remaining 
part of this work, we will limit ourselves to the determination of the static diffusion barriers. 



IV. OTHER DIFFUSION EVENTS IN THE PRESENCE OF STRAIGHT STEPS 



A. Diffusion of adatoms across straight steps: Schwoebel barrier 



When an adatom is deposited on an adisland it wanders on its surface until either it meets a group of adatoms 
already present on this adisland or it comes down the step edge by diffusing across the step. However the additional 
barrier (called the Schwoebel barrier) felt by the adatom for the latter process is often rather high relative to the 
barrier (40meV) encountered on the flat surface and the probability that an adatom crosses the step is small. As 
a consequence it is often assumed that the flux of adatoms towards the step comes from the lower terrace. The 
"traditional picture" of an adatom "jumping" down the step (process j of Fig. UJ and therefore loosing one nearest 
neighbor in the diffusion process usually leads to high energy barriers but another type of diffusion process exists, i.e., 
the exchange mechanism in which the adatom takes the place of a step atom which becomes an adatom attached to 
the step (processes xq and Xi of Fig. UJ. In Table llVl we present the energy barriers corresponding to the three types 
of diffusion mechanisms across the steps A and B on an adisland in FCC and HCP configurations. Interestingly the 
energy barriers for a jump mechanism are almost the same (around 510meV) in all cases, but the exchange mechanism 
has a lower barrier by ~ 200meV for the process xo (~ 130meV for xi) on step A and ~ 400meV for both exchange 
processes on step B. A strong anisotropy is therefore found between step A and step B. The Schwoebel barrier is thus 
only ~ 50meV for step B while it is around 250meV and 300mcV for xo and xi, respectively, on step A. The existence 
of a low Schwoebel barrier across step B is in agreement with the previous EAM calculations by Trushin et al 42j and 
could be at the origin the fast decay of double layer Cu adislands on Cu(lll) observed in experiments [3 E ■ 



B. Formation and diffusion of dimers 



A step adatom diffusing along the step will possibly meet another single adatom to form a step dimer, or stick 
to a kink. The latter case will be discussed later (see Sec.V). Let us first mention that the adsorption energy of a 
dimer along a B-step is favored by 5meV with respect to an A-step. In Table^we present the energy barriers for the 
formation and dissociation of a dimer (see FiglSJ. As expected the diffusion of an adatom is favored by the proximity 
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FIG. 7: Diffusion processes of a descending adatom across the A and B- steps iimiting a stripe in FCC geometry on Cu(lll), 
by jump (j) or exchange (xo and xi) mechanisms. The atoms are represented as in Figl2l 



Step 


A 


B 


Geometry 


FCC 


HCP 


FCC 


HCP 


Direction 


j X X! 


j X Xi 


j X X! 


j X Xl 


► 

•4 


507 299 374 
1222 1014 1032 


508 289 362 
1218 999 1016 


510 90 106 
1228 807 776 


510 86 103 
1222 799 767 



TABLE IV: Diffusion barriers of an adatom descending or ascending a step (A or B) by various mechanisms: jump (j), exchange 
(xo or xi) (see Fig. |7Jfor FCC and HCP geometries. The directions ► and -4 correspond to descending and ascending adatoms, 
respectively. 

of another adatom, the energy barrier for the dimer formation being around 15 % lower than the energy barrier to 
diffuse freely along a step. Once the dimer is formed the probability for the dimer dissociation is rather low since in 
any case the associated barrier is larger than 500meV (see Table [Vj), however a concerted motion of the two atoms is 
conceivable (see Fig . 



Step 


A 


B 


Geometry 


FCC 


HCP 


FCC 


HCP 


Dimer dissociation 


520 


509 


568 


560 


Dimer formation 


207 


197 


262 


255 



TABLE V: Energy barriers for the formation and dissociation of a step dimer along A and B-steps in FCC and HCP geometries. 

We have therefore calculated the energy barrier for the diffusion of such step dimers (Table IVlj) . Interestingly the 
barrier height for the dimer motion shows a more pronounced anisotropy between step A and step B than for the 
single adatom motion. Indeed the ratio of the diffusion barriers along steps A and B is ~ 1.27 and ~ 1.40 for the single 
adatom and concerted dimer motions, respectively. In addition the comparison of Tables IVl and IVll shows that the 
concerted motion of the dimer is easier than the dimer dissociation along step A, while along step B the dissociation 
is favored. As a consequence, the dimer diffusion will proceed by concerted motion along step A, whereas along step 
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B a two step process is preferred, i.e., a dissociation followed by the re-bonding of the dimer. 



Step 


A 


B 


Geometry 


FCC 


HCP 


FCC 


HCP 


Edim (meV) 


463 


440 


640 


625 



TABLE VI: Diffusion barriers Edim for the dimer diffusion along steps A and B in FCC and HCP geometries. 
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FIG. 8: Atomic motions corresponding to the dissociation, the formation and the concerted diffusion of a Cu dimer along the 
steps A and B limiting a stripe in FCC geometry on Cu(lll). The atoms are represented as in Fig|5| 



V. DIFFUSION OF ADATOMS ALONG STEPS WITH DEFECTS 

When an adatom sticks onto a preformed close-packed adisland bordered by A and B steps with defects (kinks 
or other adatoms) it will in general remain attached to it and, if edge diffusion and corner crossing processes are 
activated, the adatom will diffuse back and forth between the two types of steps. In the following section we will start 
by presenting the potential energy profiles for typical diffusion sequences along which the most important processes 
occur. Then the other diffusion processes will be systematically investigated and the result for their activation barrier 
will be given and commented on. 

A. Typical diffusion sequences 

The diffusion path around an adisland in FCC geometry with A and (or) B borders and the corresponding energy 
profile is shown in FigEl An adatom starting from step A passes around a kink (kink crossing), detaches from this 
kink, diffuses along step A and then passes around a corner between step A and step B (corner crossing) on which 
a similar path is followed. Such a sequence is very instructive since the most important processes are encountered 
(the formation and dissociation of dimers along the steps have been studied in Sec. IV. B): i) diffusion along step A 
(a ^ b, e ^ /)and along step B (j ' ^ k, n ^ o), ii) kink crossing on step A (6 — > c — > d) and on step B (n — > m — > I), 
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Hi) kink attachment (e — > d on step A, k — > I on step B) or detachment (d — ► e on step A. Z 
corner crossing (g ^ h ^ i). 
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FIG. 9: a) A typical diffusion path of an adatom along the borders of an adisland in FCC geometry with A and B edges on 
Cu(lll), b) The associated energy profile. The atoms are represented as in Fig|2] 

The energy barriers for the diffusion along step A and step B have been already studied (Sec. III. B). The numerical 
values given in Figl^D are very slightly different (by a few meV) due to the finite size of the adisland edges. 

The kink crossing is a two step process since a secondary energy minimum is found at site c on step A and m on 
step B. The energy barrier to pass around the kink is significantly larger than for the diffusion along straight steps. 
This is actually the one-dimensional analogue of the Schwoebel effect encountered when diffusing across steps. This 
means that the adatom is repelled by the descending kink. Furthermore even if the adatom reaches site c on step 
A, it will most probably come back to site b (barrier: 20mcV) than continue its way to site d (barrier: 85mcV). On 
the contrary on step B an adatom reaching the site m will most probably continue its way to site I (barrier: 5mcV) 
rather than turning back to site n (barrier: lllmeV). 

The energy barriers for kink detachment are rather high (505meV for step A, 559meV for step B) since roughly two 
nearest neighbor bonds have been broken at the saddle point. Moreover the energy barriers for the kink attachment 
are significantly lowered (by 30-40meV) compared with the diffusion barriers of the straight steps. Thus an adatom 
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reaching site e or k has a larger probability to attach to the ascending kink than to be reflected. It is worthwhile to 
note that this effect cannot be obtained when using approximate expressions for the barriers. Finally it is interesting 
to compare these energy barriers to those involved in the dissociation and formation of dimers along steps A and B. 
By comparing with the results of Table V, we see that the kink detachment is easier than the dimer dissociation, 
while the kink attachment is more difficult than dimer formation. This is clearly due to the change of local atomic 
environment. 

The energy profile for the corner crossing has some similarity with the kink crossing since the energy barrier is 
significantly larger than for the diffusion along straight steps, i.e., the adatom is repelled by the corner leading to a 
Schwoebel-like effect. Furthermore a secondary energy minimum is met at site h during the process. However the 
two minima on both sides of site h (sites g and i) are roughly at the same energy (the energy difference between sites 
g and i being ~ 7meV) since both sites have the same first nearest neighbor coordination. Moreover it is interesting 
to note that an adatom located at site h has a larger probability to escape towards step A (barrier: 22meV) than 
towards step B (barrier: 108mcV). 

In addition, contrary to the diffusion across the straight steps (see Table IV) the exchange mechanism is never 
favored compared with the jump process. Indeed, we have already seen (Table I) that the diffusion along steps does 
not proceed by exchange. Similarly we have verified that the exchange mechanism is also unfavorable for kink and 
corner crossings. This is illustrated in Fig. 1101 for the corner crossing: the usual double jump mechanism (g — > h — > i) 
is much preferred to a single (Fig llOb ) or double exchange (Fig llOb ) mechanism. Surprisingly however the rather 
unexpected double exchange mechanism is less costly than the single exchange. 

Bearing in mind that crystal growth is a kinetic phenomenon leading to out of equilibrium shapes like triangular 
adislands, it is also important to study typical diffusion processes occurring in the vicinity of such adisland shapes. 
For this reason we have also considered the corner crossing of a triangular adisland bordered with A or B steps. The 
results of our calculation, illustrated in Fig. 1111 is quite instructive since it is found that the corner crossing of a 
triangular adisland bordered with B steps is much easier than for a triangle with A edges. This asymmetry already 
exists if the crossing mechanism occurs by successive elementary jumps, but it is much more pronounced for the 
exchange mechanism which is highly favorable at the corner of a triangular adisland bordered with B edges. Indeed 
the corresponding activation energy is only 230meV, i.e., even much less than the simple diffusion along the straight 
step B. This could have important consequences on the growth scenario. 

B. Other diffusion events 

Finally, in view of carrying out a KMC simulation of Cu/Cu(lll) growth it is important to have a complete picture 
of the various diffusion processes. An enumeration of (almost) all diffusion mechanisms (labelled by roman figures) 
of an adatom in the vicinity of a preformed adisland is presented in Fig. 1121 For the sake of completeness we have 
also considered the possibility for a step edge atom (V, VI, VI') to escape along the step or a step adatom (IV) 
to escape on the terrace. We use the following notation hereafter: the initial configuration of the diffusion events 
represented in Fig ll2l is denoted by an index and the final configuration by an index i = 1,2,3,4. For example 
the simple diffusion of an adatom along a straight step A is denoted as: IV'^-p The arrow is simply reversed for 
the opposite displacement. Note that, even though FCC-HCP jumps are the smallest possible diffusion movements, 
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FIG. 10: Possible diffusion mechanisms for passing around the corner between step A and step B on an adisland in FCC 
geometry on Cu(lll): a)jump process, b)single exchange and c)double exchange mechanisms, d)the corresponding potential 
energy profiles: jump (full line), single exchange (dashed line) and double exchange (dotted line). The atoms are represented as 
in Fig|U 

they were not taken into consideration here (except for IIjf_^ 4 discussed in the following) since in the vicinity of an 
adisland all the sites that can be reached by such jumps are unstable, and an adatom situated on such a site would 
immediately be attracted by the step. 

Most of the diffusion processes along steps have already been discussed in detail above and we shall not make any 
further comment on them. Let us however try to extract some general trends from Table IVlTl 

One can observe that the diffusion energy of an atom in the vicinity of an adisland in HCP geometry is systematically 
smaller by a few meV than the corresponding energy in the vicinity of an adisland in FCC geometry similarly to the 
diffusion along straight steps. 

Obviously, the energy necessary for an atom to leave an adisland towards the terrace and loose all its first nearest 
neighbors from the adisland, is much higher than the diffusion energy along steps. Indeed the energy cost of such 
processes is around 950meV for a corner (IVo-^) or a kink atom (V'o->2) ; around 650meV for a step adatom to escape 
on the terrace (IV'0^2) and decreases to 400meV for the motion II^ 2 - These energy barriers obviously decrease with 
the number of bonds lost in the diffusion process which are equal to 3, 2 and 1, respectively. However this correlation 
is very approximate since, for a given number of broken bonds, the values of the activation barriers vary over a rather 
large interval: for instance, when a single bond is broken the activation barrier lies between 323 and 560meV. Indeed, 
similarly to the case of straight steps, the potential energy surface of an adatom is strongly affected by the vicinity 
of the adisland. As a consequence the diffusion path cannot be guessed from symmetry arguments and, at the saddle 
point, many different interatomic distances with the neighbors are involved so that it is difficult to derive an effective 
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FIG. 11: Upper part: possible diffusion mechanisms (jump and exchange) for passing around the corner of a triangular adisland 
bordered a) by A-steps or b) B-steps in FCC geometry on Cu(lll). The atoms are represented as in Fig|5| Lower part: the 
corresponding potential energy profiles: jump (full line) and exchange (dashed line) mechanisms. 
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FIG. 12: Diffusion processes of an adatom around an adisland bordered with A and B steps at FCC sites on Cu(lll). The 
atoms are represented as in Fial^l 

pair interaction model as done in our previous workj^lj]. Let us also stress that the terrace site that can be reached 
from an adisland site may be unstable or quasi unstable, see for example the diffusion processes IVo^2 and V'o<— 2 for 
both steps. 
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Finally let us mention that the II^_ >4 process is also affected by the presence of the adisland. Indeed the energy 
profile would be very different if the diffusing atom and its nearest neighbor were isolated on the terrace. The dimcr 
motion was studied in detail in our previous work 1 2 1 1] and the energy barrier corresponding to a motion of type 
IIo-»4 from an FCC-FCC to an FCC-HCP configuration (also called ff-fh in Refsl2ll |29|) was found equal to 16meV 
(the initial and final sites being almost energetically degenerate), whereas in the presence of an adisland this motion 
becomes highly asymmetrical, and the true dimer (FCC-FCC) is stabilized by the vicinity of the adisland. This 
phenomenon is similar with the experimental finding of Repp et al. who observed the stabilization of a dimer by 
the proximity of a monomer. In the case of II^ 4 the influence of the step atoms is so strong that the final site 
becomes unstable. 



VI. CONCLUSION 



In conclusion, we have carried out a systematic study of the diffusion processes that can occur for a Cu adatom in 
the presence of a close-packed adisland on Cu(lll). In view of the small value of the stacking fault energy in Cu, two 
geometries have been considered in which the atoms of the adisland occupy FCC or HCP sites. The diffusion rates 
along step A and step B without defects have been first investigated using TST-HA and MD simulations in the classical 
limit. They obey Arrhcnius laws with activation barriers significantly larger along step B than along step A. Whereas 
the (static) potential energy barriers account for the slope of the MD Arrhenius plots within statistical errors, the 
attempt frequencies are markedly different due to the anharmonic and recrossing effects included in MD simulations. 
Indeed, in contrast with the Vineyard attempt frequencies which have similar values for all geometries, the prefactors 
derived from MD simulations obey the Meyer-Neldel compensation rule, i.e., the increase of the activation barrier 
for B steps compared with A steps is somewhat compensated by an increase of the prefactor. Furthermore, the law 
proposed by Boisvcrt et al. jssl to relate the activation barrier to the prefactor accounts quite well for our results. 
Consequently we have then limited ourselves to the determination of static potential barriers for a large number of 
diffusion events that can occur in the presence of an adisland. A number of additional differences have been put 
forward between the diffusion along A and B borders which may have an influence on the adisland growth shapes. In 
agreement with previous EAM calculations wc find that the exchange mechanism for descending a step is largely 
favored, especially on step B for which the Schwoebel extra barrier is only ~50-60meV. Moreover it is shown that 
the corner crossing diffusion process by jump or exchange for triangular adislands with A borders have a similar and 
high (~480meV) activation barriers while for triangular adislands with B borders the exchange mechanism is largely 
favored and has a rather low barrier (230meV). From this set of results the diffusion rates of the most important 
atomic displacements can be predicted and used as input in KMC simulations which are currently in progress. 
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A 


B 


Diffusion event 




FCC 


HCP 


FCC 


HCP 


II 0-» 1 





38 


33 


189 


180 


<- 1 




38 


33 


189 


180 


0^2 


1 


400 


396 


411 


405 


0^2 




30 


28 


32 


29 


0^3 




60 


54 


3 


3 


0^3 


2 


668 


660 


600 


602 


0^4 




39 


33 


- 


- 


0^4 




6 


10 


- 
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III -> 1 


1 


327 


323 


424 


415 


<- 1 




21 


17 


113 


104 


0^2 


2 


667 


664 


689 


686 


0^2 




3 


1 


17 


16 


IV 0-> 1 


2 


609 


606 


609 


606 


<- 1 




10 


9 


10 


9 


0^2 


3 


952 


948 


952 


948 


0^2 
















0^3 


2 


695 


685 


695 


685 


0^3 




94 


87 


94 


87 


IV -> 1 





247 


235 


312 


302 


<- 1 




247 


235 


312 


302 


0^2 


2 


660 


661 


683 


677 


0^2 




0.6 


0.7 


19 


15 


V 0-> 1 


3 


913 


903 


866 


862 


<- 1 




45 


40 








V 0-> 1 


1 


506 


494 


560 


551 


<- 1 




211 


200 


270 


262 


0^2 


3 


950 


947 


947 


945 


0^2 
















0^3 


2 


681 


674 


600 


599 


0^3 




86 


80 


8 


7 


VI 0-» 1 


3 


842 


837 


824 


816 


<- 1 




45 


43 


19 


16 


vr o i 


2 


726 


713 


647 


633 


<- 1 




196 


183 


114 


103 



TABLE VII: Calculated diffusion barriers for a Cu adatom jumping from an initial site to a final site (1, 2, 3, or 4) and 
vice-versa. AZ is the number of nearest neighbor bonds broken in the diffusion process. The various diffusion events are shown 
in Fig[ni 
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APPENDIX A 

Let us consider the (3N+3-dimensional) configuration space of a system with N+f identical atoms of mass m and 
assume that the potential energy of this system has two minima at points A and B separated by a saddle point P. We 
suppose that the dividing surface £p (containing P) between the regions A and B as well as the diffusion path (i.e., 
the steepest descent line going from A to B through P) have been determined (see Fig li3(l . The coordinate system is 
chosen so that one coordinate, s, runs along the diffusion path. The equation of Ep is then s = sp. In the transition 
state theory (TST) the diffusion rate from A to B is given by Q: 

t^A — >B /rji\ _ < V s 0{Vs)S(s - S P ) > 

1tst(J) ~ <9(s p -s)> (A1) 

where v s is the velocity along s, < ... > denotes a canonical average and 9 is the Heaviside function. 

The basic assumption leading to Ea lAll is that each crossing of Sp corresponds to a diffusion event. This is not 
true since recrossing may occur before thermalization in region B. To correct for this effect a transmission coefficient 
k{T) is introduced so that the total diffusion rate is: 



r(T) = n c K{T)T^*(T) (A2) 

when there are n c diffusion channels. 

1. Diffusion rate in the TST-HA 

When integrating over p s = mv s in the canonical average of the numerator of Ea lAf I T^^p (T) can be rewritten: 

rrST B (T) = ^f^ (A3) 

where Zs p is a constrained partition function, i.e., in which the representative points of the system in the 3N+3- 
dimensional configuration space arc compelled to stay on Ep and Z a is the partition function when the representative 
points are located in the region around A. Thus the dimensionality of the configuration space is 3N+2 for Z^ p and 
3N+3 for Z A . 

If the static barrier AE is much larger than k B T, the HA can be used and r^^?(T) becomes: 

^tst(T) = C M-(&E + AF mb )/k B T) (A4) 
where AF v n, = AU v a, — TASmt is the contribution of vibrations to the free energy, with: 

a tt r v maa fa fa 

AU lnb = / — coth (— — — )An(v)dv 
Jo 2 x 2k B T' 



ASyib — fcp 







cothf ) - In (2s nh( )) 

2k B T K 2k B T> V y 2k B T JJ 



An(v)dv (A5) 
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and: 



An(v)du = -1. (A6) 

Indeed, in the HA the constraint on Z^ p can be expressed as qp = where qp is the normal coordinate corre- 
sponding to the unstable mode at point P. As a consequence this mode must be excluded from the summation over 
v. Consequently, in this approximation and in the classical limit: 



AS v ib 



k B T 

Kb 2 



(A7) 
(A8) 



where Vq is the Vineyard attempt frequency fEa lll|) . These quantities have been calculated explicitly for self-diffusion 
on the (100) surface of Cu and Ag |3lJ,|32|. Finally note that the correction factor k(T) is generally omitted in the 
HA. 

2. Diffusion rate from the TI method 

In this approach, the integration over all momenta in both canonical averages of Ea lAll is carried out in the classical 
limit. This leads to 



1 27T771 ' 



f exp{-E{...n j ...)/k B T)dZp 

1/2 JSp 

f exp(-E(...r tJ ...)/k B T)dV A 
Jv A 



(A9) 




FIG. 13: Schematic potential energy surfaces in the vicinity of two potential wells A and B separated by a saddle point P. The 
surfaces E s are given by s=cst where s is the coordinate along the diffusion path. 

The integration domain extends over the hyper-surface Ep in the numerator and the hyper-volume Va limited by Ep 
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and containing A in the denominator. Let H s be the surface s = est, i.e., it is normal at s to the steepest descent line 
containing A and P and set: 



Cexp(-W(s)/k B T) 



cxp(-E(...r lj ...)/k B T)dT, s 



(A10) 



where C is a constant. If we multiply the numerator and denominator of Eci lA9l bv exp(— W{sA)/k B T), this equation 
can be rewritten as: 



27T7T1 

= v{T) exp(- AW/k B T) 



P exp(-(W{a) - W{s A ))/k B T)ds\ ' exp(-(W(s P ) - W{s A ))/k B T) 

— OO 



(All) 
(A12) 



The constant C can be fixed by setting W(s A ) = 0. Thus 



W{s) = -k B T\n 



J exp{ r E(...r ij ...)/k B T)dS t 
[ exp{-E(...r ij ...)/k B T)dZ J 



(A13) 



The function W{s) is known as the "potential of mean force" 56] in the literature. Indeed it is straightforward to 
show that 



dW 
ds 



dE 
ds 



(A14) 



Finally AW = W(sp) — W{sa) is an activation free energy which satisfies: 



exp{-AW/k B T) = 



(A15) 



in which both Z^ p and Z^ A are now constrained partition functions having thus the same dimensionality (3N+2) of 
the configuration space. Note that this point of view is closely related to the formulation of Wert and Zener |47l l57j. 
We must emphasize that if we set AW = AE + AU TI - TAS TI , the quantities AU TI and AS TI are different from 
AU V ib and AS V ib (Ea lA8|l even in the HA. Indeed the number of modes, being the same in the two partition functions, 
AUjj A vanishes and AS]^ is independent of temperature. 

This method has been used by Boisvert et a/.|3jj( for the diffusion of Cu on Cu(100) without resorting to the HA. 
They used MD simulations to calculate AW, v(T) and k(T) and showed that their result for T(T) fEa lA2(l can be 
fitted nicely by an Arrhenius law. Indeed they found that: i) if AW is written as AE TI — TAS TI , then AE TI and 
AS TI are both effectively temperature independent, ii) k(T) and v(T) are only slightly dependent on temperature. 
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